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Dear  Jon, 

As  a  result  of  several  conversations  including  one  with  Richard  Orr 
of  Atlantic  Aerospace  I  have  become  concerned  with  the  efforts  of 
Ana  Tsao,  Lou  Auslander  and  others  propagandizing  for  the  Auslander- 
Johnson  FFT  algorithm  to  exagerate  their  efforts  at  the  expense  of 
mine  and  those  of  my  students  and  collaborators.  I  am  simply 
interested  in  protecting  the  priority  of  these  efforts  and  in  no 
way  commenting  on  the  Auslander-Johnson  efforts  as  they  involve < 
the  integrity  and  influence  peddling  of  DARPA  contracts  over  which 
I  have  no  control. 

Inclosed  are  two  papers:  One  published  by  my  group  in  1991  and  the 
second  recently  submitted  by  Auslander  et  al  in  1996  which  proclaims 
the  Auslander-Johnson  FFT  supported  by  DARPA. .  Both  describe  the 
exact  same  divide-and-conguer  algorithm  for  finite  abelian  groups. 

The  relevant  pages  highlighted  in  yellow  are  pages  25-26  in  my  work 
and  pages  7  and  8  in  Auslander 's  work.  Note  the  Auslander  work  is 
equivalent  to  his  previous  efforts  involving  induced  representations 
as  stated  on  page  3.  Significant  extensions  of  this  approach  have 
been  completely  described  in  my  second  Springer  book  and  in  a  second 
paper  also  enclosed  relating  the  abelian  group  FFT  to  affine  group 
actions. 

These  works  have  been  submitted  along  with  implementations  to  both 
AFOSR  and  DARPA  from  1989,  and  have  been  mailed  and  hand  delivered 
to  both  Auslander  and  Tsao.  In  addition,  tensor  product  construction 
for  the  FFT  matrix  multiplication  and  the  DCT  along  with 
implementations  and  automatic  code  generation  have  been  carried  out 
in  1989  by  two  of  my  students  John  Granata  and  Martin  Rof heart 
who  currently  run  their  own  businesses  in  DC  aera. 

Note  also  my  paper  is  a  tutorial  survey  which  only  claims  to 
mathematically  organize  the  works  of  others  which  makes  it  even 
more  ironic  that  the  new  Auslander-Johnson  FFT  algorithm  can  be 
thought  of  as  new  five  years  after  publication  of  the  survey. 
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Review  of  "Radar  Waveform  Design  and  Clutter  Suppression"  by  L.  Auslander 

Dr.  Auslander  proposes  to  design  radar  waveforms  to  suppress  clutter  in  SAR  images.  However,  he 
does  not  explain  why  he  wants  to  do  this!  It  is  not  clear  whether  he  really  wants  to  suppress  clutter  or 
sidelobes.  If  sidelobe  suppression  is  the  goal,  then  the  SVA  algorithm  in  Reference  [1]  is  the  best 
method  (as  far  as  I  know),  but  it  is  not  mentioned  in  Dr.  Auslander's  proposal. 

The  first  two  sentences  of  the  section  titled  "Scientific  Proposal"  on  Page  1  state:  "Classical  clutter 
suppression  at  the  radar  waveform  level  -  say  for  Doppler  -  has  been  achieved  by  studying  the 
spectrum  of  the  clutter  and  designing  waveforms  whose  Fourier  transforms  have  little  energy  in  the 
support  of  the  spectrum  of  the  clutter.  Classical  clutter  suppression  at  the  level  of  SAR  images  has 
been  based  on  statistical  studies  of  the  clutter  and  the  design  of  pre-whitening  filters." 

"Clutter"  is  clearly  undesirable  in  MTI  radar  applications.  However,  in  SAR  applications  the  clutter  is 
the  signal!  Saying  that  you  want  to  suppress  the  clutter  in  SAR  images  is  equivalent  to  saying  that 
you  want  to  suppress  the  signal  of  interest,  unless  you  carefully  explain  what  you  mean  by  clutter 
suppression.  Unfortunately,  I  could  not  find  this  explanation  in  the  proposal/ 

Clutter  suppression  is  usually  considered  to  be  a  disadvantage,  not  an  advantage,  in  SAR  super¬ 
resolution  algorithms.  These  algorithms  achieve  resolutions  better  than  (1/Bandwidth)  at  the  expense 
of  emphasizing  isolated  point-like  discretes  over  the  more  complicated  clutter.  The  details  of  diffuse 
clutter  features,  such  as  vehicle  tracks,  are  very  important  in  SAR  reconnaissance  applications. 
Algorithms  that  suppress  these  clutter  features  are  usually  considered  to  be  bad,  not  good! 

There  may  be  some  benefits  to  suppressing  diffuse  clutter  in  imagery  that  will  be  analyzed  by  an 
automatic  target  recognition  (ATR)  algorithm,  but  this  application  is  not  mentioned  in  Dr.  Auslander's 
proposal.  (If  Dr.  Auslander  is  indeed  tTying  to  suppress  diffuse  clutter  for  an  ATR  algorithm,  then  he 
would  need  to  be  specific  about  which  ATR  algorithm  he  plans  to  use.) 

Dr.  Auslander's  mathematics  may  lead  to  exotic  waveforms  that  practical  hardware  cannot  generate  or 
propagate.  Very  high-resolution  SAR  images  requiring  very  wideband  waveforms  are  usually 
obtained  with  linear  FM  waveforms  because  it  is  difficult  to  design  receiver/exciters  and  high-gain 
directive  antennas  with  other  types  of  very  wideband  waveforms.  This  practical  issue  is  not 
addressed  in  the  proposal. 

The  example  provided  in  the  proposal  is  based  on  ECG  data  for  heart  beats.  Radar  examples  would 
have  been  much  more  appropriate. 

In  summary,  my  impression  of  the  proposal  is  very  negative. 

Reference 

[1]  H.C.  Stankwitz,  R.J.  Daillaire,  J.R.  Fienup,  "Nonlinear  Apodization  for  Sidelobe  Control  in 
SAR  Imagery,"  IEEE  Transactions  Aerospace  and  Electronic  Systems,  January  1995. 
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Abstract1 —  This  paper  presents  a  divide  and  conquer  algorithm  to  compute 
the  Fourier  transform  of  a  finite  Abelian  group  using  the  Fourier  transform  of  an 
arbitrary  subgroup  and  the  Fourier  transform  of  the  corresponding  quotient  group. 
The  construction  uses  an  arbitrary  choice  of  coset  representatives  for  the  quotient 
group.  Different  choices  of  coset  representatives  lead  to  different  data  flows  in  the 
algorithm  and  different  twiddle  factors. 

The  derivation  of  the  algorithm  generalizes  the  derivation  of  the  FFT  presented 
by  Cooley  and  Tukey.  Moreover,  it  can  be  used  to  obtain  explicit  algorithms  for 
multidimensional  Fourier  transforms.  The  algorithm  presented  in  this  paper  gener¬ 
alizes  all  of  the  known  Cooley-Tukey  type  algorithms  for  multidimensional  Fourier 
transforms  and  provides  new  algorithms  with  alternative  data  flow  patterns.  Some 
preliminary  experiments  suggest  that  in  hierarchical  memory  computers  these  al¬ 
gorithms  are  more  efficient  than  the  standard  “row-column”  approach. 
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1  Introduction 

In  1965,  Cooley  and  Tukey  [11]  presented  a  divide  and  conquer  algorithm, 
called  the  Fast  Fourier  Transform  (FFT),  for  computing  the  Fourier  trans¬ 
form  of  a  sequence  of  n  points.  The  FFT  has  a  long  and  interesting  his¬ 
tory  [10]. 

Cooley  and  Tukey  considered  the  problem  of  computing 

X(j)  =  £  ;  =  0,1, -  1,  W  =  e7**tN  (1) 

k=Q 

given  a  complex  function  A(k). 

Their  algorithm  is  obtained  by  viewing  the  functions  X  and  A  as  functions 
of  two  variables  and  rearranging  the  sum  in  Equation  1  as  a  nested  sum.  Since 
the  material  in  this  paper  is  closely  modeled  on  the  ideas  in  the  paper  by 
Cooley  and  Tukey,  we  reproduce  their  derivation. 

To  derive  the  algorithm,  suppose  N  is  composite,  i.e.,  N  —  rir2.  Then  let 
the  indices  in  (1)  be  expressed 

j  =  /[Ti  +  Jo  I  jo  =  0,  1  ■  •  •  •  >  Cj  —  1,  j\  —  0,  1,  •  .  -  I  *2  —  1|  0^) 

Ic  —  Tn  *f  ^>‘o  *  ^0  —  0|  1 1  •  ■  •  (  ^1  —  0,  1|  •  *  •  I  rl 

Then,  one  can  write 

A'O'i.jo)  =  (3) 

ko  k  | 

Since  .  x 

\yjkxr7  _  \y)okxri'  (4) 

the  inner  sum,  over  k\ ,  depends  only  on  jo  and  to  and  can  be  defined  as  a 
new  arrav, 

Al(j0.k0)  =  Y,A^'k^wittir3-  (5) 

kx 

The  new  result  can  then  be  written 

A'Uiiio)  =  Y,AMMwU,ri+io)ko-  (6) 

ko 

Alternative  algorithms  have  been  presented  when  ri  and  t-i  are  relatively 
prime  [13].  The  approach  of  Cooley  and  Tukey  has  been  used  to  obtain  mul¬ 
tidimensional  Fourier  Transform  Algorithms.  The  “vector  radix”  algorithm 
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was  presented  in  [15,  23],  and  a  more  general  multidimensional  algorithm 
was  independently  presented  in  [4]  and  [22]. 

Recently  the  authors  have  generalized  the  Cooley-Tukey  algorithm  to 
apply  to  an  arbitrary  finite  Abelian  group  [1].  This  algorithm  computes  the 
Fourier  Transform  of  a  finite  Abelian  group  using  the  Fourier  transform  of  an 
arbitrary  subgroup  and  the  Fourier  transform  of  the  quotient  group  obtained 
by  moding  out  by  the  chosen  subgroup.  The  construction  uses  an  arbitrary 
choice  of  coset  representatives  for  the  quotient  group.  Different  choices  of 
coset  representatives  lead  to  different  data  flow  in  the  algorithm  and  different 
twiddle  factors.  These  results  generalize  all  of  the  known  Cooley-Tukey  type 
algorithms,  including  those  mentioned  above,  and  provides  new  algorithms 
with  alternative  data  flow. 

The  results  in  [1]  were  obtained  using  the  theory  of  induced  representa¬ 
tions;  similar  in  spirit  to  recent  work  on  algorithms  for  the  computation  of 
Fourier  transforms  of  non-Abelian  groups  [S,  9,  12,  24].  In  this  paper  an  al¬ 
ternative  combinatorial  proof  is  presented.  This  proof  generalizes  the  original 
proof  used  by  Cooley  and  Tukey.  The  benefit  of  this  approach  is  that  it  pro¬ 
vides  an  explicit  formula  that  can  be  used  for  constructing  multidimensional 
Fourier  Transform  algorithms.  A  meta-algorithm  is  provided  that,  given  a 
presentation  for  an  Abelian  group,  a  presentation  for  a  subgroup,  a  set  of 
coset  representatives  for  the  quotient  group,  and  a  set  of  coset  representatives 
for  a  related  quotient  group  in  the  dual  group,  constructs  a  matrix  factor¬ 
ization  of  the  Fourier  transform  of  A.  This  matrix  factorization  provides  an 
algorithm  for  computing  the  given  multidimensional  Fourier  transform. 

In  Section  2  the  Fourier  transform  of  a  finite  Abelian  group  is  denned.  A 
proof  of  the  main  result  is  presented  in  Section  3.  Section  4  shows  how  to 
apply  the  results  about  Abelian  groups  to  multidimensional  Fourier  trans¬ 
forms,  and  Section  5  presents  some  concrete  examples.  In  particular,  it  is 
shown  how  to  obtain  the  existing  algorithms  as  special  cases  of  the  gen¬ 
eral  algorithm.  In  addition  a  class  of  examples  are  presented  that  could  not 
have  been  computed  with  existing  algorithms.  Finally,  Section  6  discusses 
the  consequences  of  the  results  in  this  paper  to  the  implementation  of  high 
performance  multi-dimensional  Fourier  transform  algorithms. 
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The  Fourier  Transform  of  a  Finite  Abelian 
Group 


Let  A  be  a  finite  Abelian  group  whose  order  is  denoted  by  \A\  (the  group 
operation  will  be  denoted  by  +).  The  set  of  complex  valued  functions  on  A 
with  inner  product 

(/,$)  =  777 

\A\  ae A 

form  an  inner  product  space  denoted  by  L2(A). 

A  non-zero  function  x  €  L2(A)  such  that 


x(ai  +  c2)  =  x(ai)x(a2)  f°r  ®i> ®2  €  A, 

is  called  a  character.  If  aJz_A  is  of  order  n,  then  x(°)n  =  !>  and  therefore 

lx(o)l  =  1  and  x(-a)  =  x(a)-  .  . 

If  xi  and  X2  are  characters  then  we  can  define  (xi  +  X2)(a)  —  Xi(aM2Vaj- 
Under  this  operation  the  set  of  characters  form  a  group,  called  the  dual  or 
character  group,  which  we  will  denote  by  A.  It  is  well  known  that  A  is 
isomorphic  to  A  (see  Section  4  and  [14]). 

Suppose  that  x  ^  1,  then  there  exists  Gi  €  A  such  that  x(®i)  5^  b  Since, 

*(a)  =  11,  x(a  +  =  x(°)> 

agA  agA  °€A 

(1  -x(ai))  EagA  x(a)  =  0,  ?-nd  since  x(ai)  #  1.  EogA  x(a)  =  0-  The  following 
property  easily  follows  from  this  calculation. 


Lemma  2.1  Let  Xi:X2  €  A. 

(X11X2) 


f  1  if  Xi  =  X2 
\  0  if  Xi  X2 


This  lemma  implies  that  A  is  an  orthonormal  basis  for  L2(A).  Hence,  an 
arbitrary  function,  /  €  L2{A)  can  be  written  uniquely  as 

/(*)  =  J2  /(x)x(*).  where  fix)  =  (/.x)-  (7) 

x€A 

This  is  the  Fourier  series  expansion  of  the  function  /.  The  coefficients  in 
this  expansion  are  called  the  Fourier  coefficients  and  are  obtained  from  the 
Fourier  transform  of  /. 
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Definition  2.1  (Fourier  Transform) 

The  Fourier  transform  of  A,  F(A)  :  L2(A)  —  L7{A)  is  defined  by 

F(A)(f)(X)  =  f(x)  =  ^j£  /(fl)xP)  =  (/•  X)- 

3  A  Divide  and  Conquer  Algorithm 

Let  ( ,  )  :  A  x  A  -»  C  be  the  bilinear  pairing  of  A  with  its  dual  A  defined  by 

(a,a)  =  a(a). 

Let  B  <  A  be  a  subgroup  of  A.  Corresponding  to  B  we  can  form  a  subgroup, 
of  A  equal  to  {a  6  A\  ( b,a )  =  l/for  b  €  B)t  the  characters  that  are 
perpendicular,  with  respect  to  (, ),  to  the  subgroup  B. 

Since  a  €  B1  is  equal  to  1  on  B,  we  can  identify  it  with  a  character  on 
C  =  A/B,  by  setting  a(a  +  5)  =  a(a).  Likewise  given  a  character,  c,  of  A/B, 
we  can  obtain  a  character  in  B1,  by  composing  the  projection  A  -»  A/B 
with  c.  Therefore  B 1  is  isomorphic  to  C. 

A  character,  a  €  A,  restricted  to  the  elements  of  B  is  a  character  of  B 
denoted  by  a\B.  The  restriction  map,  a  ■->  a\B  is  a  homomorphism  from 
A  _+  B  with  kernel  B x.  Since  (B1!  =  \A/B\  and  |A|  =  |A|,  \A/ B^\  =  |B| 
and  the  restriction  map  is  onto;  therefore 

A/B1  =  B. 

As  a  result  of  this  isomorphism  there  are  I#1!  =  |C|  characters  in  A  that 
restrict  to  each  character  in  B.  Moreover,  these  are  the  characters  in  the 
cosets  of  A/B1. 

Let  C  =  A/B  and  B  =  A/B1.  Let  (  :  C  -*  A  and  rj  :  B  -»  A  be  a  choices 
of  coset  representatives.  The  following  theorem  shows  how  to  compute  F{A) 
using  \C\  copies  of  F(B),  |B|  copies  of  F{C)  and \A\  complex  multiplications. 
The  complex  multiplications  depend  on  the  choices  of  coset  representatives 
£  and  rj. 

To  simplify  the  notation  used  in  the  following  theorem  and  proof  we  will 
remove  the  normalization  constant  1/|A|  and  conjugation  in  the  definition  of 
the  Fourier  transform  and  compute 

F(A)(d)  =  £/(a)<a,a). 
oe.4 
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There  is  no  harm  in  doing  this  since  the  constants  and  conjugation  can  easily 
be  reinserted. 

Theorem  3.1  Let  B  <  A,  C  =  A/B,  C  =  Bx,  B  =  A/C,  (:  C  —*  A  be  a 
choice  of  coset  representatives  for  A/B,  and  rj :  B  —*  A  be  a  choice  of  coset 
representatives  for  A/C.  Then 

/(a)  =  £(c’£)  f(£(cM(^))  £/«<#)(*•*))  ’ 

e  €C  V  *€£>  / 

where  f^c^[b)  =  f(b  +  £(c))  and  a  €  A  ~  77(6)  +  c,  with  b  G  B  and  c  G  C. 

Proof.  Using  the  coset  decompositions  of  A/B  and  A/C  the  indexing  sets 
can  be  written  as  A  —  B  x  £(C)  and  A  =  C  x  and 

/(«)  =  ]C/(a)(a>«) 

=  EIfa(W  +  fW.H*i».  (9) 

cecteB 

Using  the  bilinearity  of  <,>  this  is  equal  to 

E  E  /«.>(*)  <*.c)  «<c),  c)  (6,4(i))  m.m).  do) 

C6C6€S 

Since  c  €  f?x,  <  6, c  >=  1  and  Equation  10  is  equal  to 

E  E  (U) 

cecteB 

Furthermore,  since  (£(c),c)  and  (£(c),t/(6))  do  not  depend  on  the  inner  sum¬ 
mation  index  6,  this  is  equal  to  the  nested  sum 

H<e(c),c)  f(€(c),i}(6))  E  /««,W(M(S)> )  •  (12) 

c6C  \  6gB  / 

Finally,  since  (£(c).c)  and  (6,77(6))  are  independent  of  the  choice  of  coset 
representatives  £(c)  and  7/(6)  this  completes  the  proof.  I 

This  theorem  is  the  basis  for  a  divide  and  conquer  algorithm  for  comput¬ 
ing  F(i4)1  which  we  now  describe. 


6 


1.  Compute  /{(c)  =  F(B)fnc)  for  c  e  C. 

2.  Compute  <7-(i)  for  b  6  B,  where  ^(c)  =  (((c),  j ?(fc))/{(c)(&)  for  c  €  C. 

3.  Compute  g.(i)  =  F(C)^(i)  for  6  €  B. 

Observe  that  /^(c)  €  L2(B),  g^  €  L2(C),  g ^  €  L2(C),  and  by  Theorem  3.1 
g.(i)(c)  =  f{fi(b)  +  c).  Since  any  a  can  be  written  as  17(6)  +  c  for  some  be  B 

and  c  e  C,  all  of  the  values  of  the  function  /  =  F{A)f  have  been  computed. 

To  implement  this  algorithm  the  elements  of  A  and  A  must  be  ordered. 
Ordering  the  elements  of  A  and  A  fixes  the  representations  of  the  functions  / 
and  /,  introducing  data  permutations  or  index  computations  when  accessing 
the  values  of  /€(c)  and  g.[iy  and  storing  the  values  of  <7^. 

If  the  elements  of  A  are  ordered  then  /  can  be  viewed  as  a  vector  whose 
indices  are  the  elements  of  A.  Step  (1)  of  the  above  algorithm  first  creates 
|C|  functions  of  B  denoted  by  /{(c).  If  the  elements  of  B  are  ordered  then 
the  functions  f^c)  can  be  represented  by  subvectors  obtained  from  /.  The 
subvectors  f^c)  can  be  ordered  using  the  order  of  the  coset  representatives 
((C).  Therefore  the  first  part  of  Step  (1)  corresponds  to  a  permutation  of 
the  vector  /.  This  permutation  is  determined  by  the  order  of  the  elements 
of  A ,  the  subgroup  B,  and  the  coset  representatives  ((C).  It  corresponds  to 
permuting  the  indexing  set  A  to  ((C)  x  B. 

After  the  vectors  f^c)  are  gathered,  f^c)  =  F(B)f^(c)  is  computed  for 

each  c  €  C.  Step  (2)  then  creates  |B|  =  |B|  functions  of  C.  For  each  be  B 
the  function  defines  a  function  of  C.  If  the  functions  f(,(c)  are  stored  as 
vectors  the  computation  of  g^i)  requires  a  stride  permutation  followed  by  a 
diagonal  multiplication  (called  the  “twiddle  factor”).  The  stride  permutation 
gathers  the  elements  of  each  vector  /€(c>  indexed  by  b.  The  twiddle  factor 

multiplies  ft(c){b)  by  {Z{c),tj(b)). 

Finally,  Step  (3)  performs  \B\  computations  of  the  Fourier  transform  on 
C,  g.(-b)  =  F(C)g^iy  As  mentioned  above,  the  resulting  functions  combine 

to  give  all  of  the  values  of  /;  however,  the  values  are  not  necessarily  in  the 
order  specified  for  A.  Instead,  the  functions  are  indexed  by  17(B)  x  C,  and  a 
permutation  similar  to  the  one  in  Step  (1)  is  required  so  that  the  resulting 
function  is  ordered  corresponding  to  the  order  of  A. 
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Abstract 

For  functions  with  finite  time  and  frequency  energy  moments,  we  find 
upper  bounds  for  the  error  of  finite  Fourier  transform  approximation  to  the 
Fourier  transform.  The  error  can  be  measured  as  the  maximum  error  over 
all  of  the  points  of  the  FFT,  or  using  a  discrete  L3  distance,  or  using  a 
continuous  L3  distance. 

Using  the  machinery  developed,  we  also  find  an  upper  bound  for  the 
error  of  approximating  the  L3  norm  of  a  function  by  a  Riemann  sura.  From 
this  result,  an  upper  bound  is  also  derived  for  the  error  of  approximating 
the  L3  inner  product,  as  well  as  the  error  of  approximating  the  integral  of 
an  L1  function,  by  a  Riemann  sum. 

As  an  application  of  the  L3  norm  approximation  theorem,  we  prove 
an  analog  of  the  Landau-Pollak-  Slepian  approximate  dimension  theorems 
for  a  certain  set  of  functions  that  is  approximately  time-and-bandlimited 
for  large  duration  N  and  bandwidth  M.  This  set  can  be  approximately 
parameterized  with  NM  parameters,  with  the  error  approaching  zero  as 
NM  approaches  infinity. 
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0  Introduction 

One  of  the  most  important  operations  in  the  signal  processing  industry  is  taking 
the  Fourier  transform  of  a  function.  With  digital  signal  processing,  the  finite 
Fourier  transform  (FFT)  is  used  instead.  A  question  of  great  practical  importance 
immediately  arises:  How  good  an  approximation  is  the  FFT? 

The  Fourier  transform  at  a  particular  frequency  is  nothing  more  than  an  in¬ 
tegral  And,  the  FFT  at  a  particular  frequency  is  a  special  case  of  numerical 
integration  using  Riemann  sum  approximation.  It  would  not  be  surprising  if  the 
theory  of  Riemann  sum  approximation  allows  us  to  derive  an  error  bound  for  the 
FFT.  In  fact,  we  will  do  just  the  opposite.  We  will  use  an  error  bound  for  the 
FFT  to  derive  an  error  bound  for  more  general  Riemann  sum  approximation. 

To  obtain  any  results,  it  is  clear  that  the  class  of  functions  must  be  restricted 
somehow.  For  example,  if  there  are  no  restrictions  on  the  oscillations,  knowing 
the  function  at  a  finite  set  of  points  tells  us  nothing  about  the  function  anywhere 
else.  We  will  consider  functions  that  decay  sufficiently  fast  in  time  and  frequency 
(in  the  sense  that  will  be  described  below).  For  this  class  of  functions,  we  will 
derive  an  error  bound  for  FFT  approximation  to  the  Fourier  transform  and  an 
error  bound  for  more  general  Riemann  approximation. 

0.1  Background 

Previous  Results 

All  previous  error  bounds  that  I  am  aware  of  come  out  of  the  observation  that 
we  can  rearrange  the  terms  of  the  Poisson  summation  formula  and  the  formula 
becomes  the  desired  error  bound.  It  is  not  clear  who  first  discovered  this.  See 
Butzer-Stens  [BS]  and  Briggs-Henson  [BH,  Chapter  6]. 

Assume  that  /  satisfies  whatever  conditions  are  necessary  so  that  the  Pois¬ 
son  summation  formula  is  valid.  Consider  the  case  of  approximating  the  Fourier 
transform,  denoted  /fr),  by  the  FFT  at  7  =  0.  That  is,  we  will  approximate 

.  ,  NMf  2-1 

/(«)  =  /  fm  by  -  £  f(k/M). 

M  k=-SM/2 

N  and  M  will  be  assumed  throughout  this  introduction  to  be  even  positive  num¬ 
bers.  N  is  the  timelength  of  the  approximation  and  M  is  the  sampling  rate. 

Note  1:  By  focusing  on  what  seems  to  be  just  the  special  case  of  FFT  approxi- 
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mation  where  7  =  0,  we  are  actually  considering  the  general  case  of  Riemann  sum 
approximation. 

Note  2:  For  the  case  of  7  ^  0,  we  can  apply  the  following  method  to  ft~7rxtl. 
In  fact,  by  the  same  method,  we  can  obtain  an  upper  bound  for  the  maximum 
error  over  all  the  points  upon  which  the  FFT  is  defined.  (See  Section  5.1.) 

One  form  of  the  Poisson  summation  formula  is 

=  i  £/{*/«)• 

z  Jeez 

A  rearrangement  of  the  terms  yields  the  following  error  bound: 

/(o) -5  E  /(*/«)  <  E  +  £  )/(*/»)  + 

Jfes-lNAf  k<-\NM  k>\NM-\ 

The  error  bound  can  be  interpreted  as  follows.  Discretization  results  in  two 
types  of  error.  Truncation  error  comes  from  the  fact  that  the  FFT  only  utilizes 
function  values  over  a  finite  duration.  Aliasing  error  comes  from  the  fact  that 
the  FFT  is  limited  to  a  finite  sampling  rate.  The  first  term  of  the  error  bound  is 
primarily  caused  by  the  truncation  error.  The  second  term  is  primarily  caused  by 
the  aliasing  error. 

The  form  of  the  error  bound  provides  practical  information  about  how  the 
error  can  be  reduced.  Let’s  say  the  first  term  of  the  error  bound  is  much  greater 
than  the  second.  Then  increasing  N  will  be  much  more  effective  in  reducing  the 
error  than  increasing  M. 

At  this  point,  we  can  impose  various  conditions  on  /  to  obtain  many  different 
theorems.  One  possibility  is  as  follows.  Let  j  be  an  integer  greater  than  1 .  Let’s 
assume  that  /k-1>  is  absolutely  continuous  (which  implies  that  fW  exists  almost 
everywhere)  and  that  /!*!  €  L'(R)  for  all  £  <  j. 

Then,  using  integration  by  parts  j  times. 

Hi) = / 

We  will  show  in  the  next  paragraph  that  the  boundary  terms  from  all  of  the 
integrations  by  parts  are  zero. 

Consider  the  following  two  typical  terms  involved  in  one  of  the  the  above  j 
integrations  by  parts: 

lim  {^)k  [B  fW(t)e-79^dt  and  lira  (~-)k+1  /*  /(fc+l)(0c'3'*^- 
il^-o.'2«7/  J-A  A,B—eo  2irif  J-A 


£  Km 
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Because  of  our  L1  assumption,  these  limits  exist  for  any  fixed  7  and  all  0  <  k  < 
j  —  1.  Therefore,  the  boundary  term  from  each  of  the  integrations  by  parts  must 
also  approach  a  limit  as  At  B  approaches  infinity.  The  boundary  term  is 

Jfe.  ("<2 ■ 

Because  of  the  exponential  factors,  the  only  way  that  this  limit  can  exist  is  if 
lim|t|— oo  f^(i)  —  0.  Therefore,  the  boundary  term  also  equals  zero. 

From  the  expression  for  /  in  terms  of  we  see  that  /  is  0(|7|-J).  Now,  let’s 
estimate  the  second  error  term  in  the  rearrangement  of  the  Poisson  summation 
formula.  The  decay  of  /  tells  us  that  the  second  error  term  is  bounded  by  a 
constant  times  £(fcM)~'7,  which  is 

We  can  also  make  the  same  assumptions  for  /  that  we  did  for  /.  Assume  that 
/D-i)  is  absolutely  continuous  and  that  6  L*(R)  for  all  k  <  By  similar 
reasoning,  we  see  that  /  is  0(|f|-;).  Therefore,  the  first  error  term  is  bounded  by 
a  constant  times  /j which  is 


L3  Approach 

The  prior  results  show  that  if  the  derivatives  of  /  and  /  are  absolutely  continuous 
and  in  L1(R),  then  we  can  obtain  error  bounds  for  both  FFT  approximation  and 
Riemann  sum  approximation.  What  happens  if  we  replace  the  hypothesis  that 
the  derivatives  are  in  L1  with  the  hypothesis  that  the  derivatives  are  in  L2?  That 
is  the  question  that  we  will  address. 

Note  that  we  are  weakening  the  hypothesis  (on  the  first  j  —  1  derivatives). 
We  saw  above  that  if  /W  is  absolutely  continuous  and  /M, /(*+*)  g  L1,  then 
limjti-,*  =  0*  This  implies  that  |/^(t)|2  <  |/^(0I  for  sufficiently  large  t 

which  implies  that  /W  €  L2.  Since  our  hypothesis  is  weaker,  we  should  expect 
weaker  results. 

There  is  a  simpler  way  to  describe  the  conditions  in  our  hypothesis.  The 
assumption  that  fb~l)  is  absolutely  continuous  and  that  /M  G  L2(R)  for  all 
k  <  j  is  equivalent  to  the  assumption  that  the  jth  frequency  moment  is  finite. 
This  moment  is  defined  by 

D)  =  /  Y'l/1^7. 

This  result  is  essentially  the  content  of  Lemmas  1.3  and  1.4. 

The  dual  of  this  result  is  also  true.  The  assumption  that  /k-1)  is  absolutely 
continuous  and  that  /W  6  L2(R)  for  all  k  <  j  is  equivalent  to  the  assumption 
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that  the  jtk  time  moment  is  finite.  This  moment  is  defined  by 

c]  =  J  t’Wdt. 

This  equivalence  may  make  our  results  more  practical  than  the  results  obtained 
with  the  L1  hypothesis.  With  the  Ll  hypothesis,  the  error  bound  depends  on  the 
L1  norms  of  the  jtk  derivatives.  This  information  is  usually  difficult  to  calculate. 
Our  approach  uses  energy  moments,  which  may  be  easier  to  calculate. 


0.2  Notation  and  Definitions 


The  version  of  the  Fourier  transform  that  we  use  is 

/b)  =  Jim  (in  V)  t 

J-A 

The  limit  refers  to  convergence  in  L2(R). 

||/||  denotes  the  norm  of  /  in  L2(R).  The  La([— N/2,  N/2])  norm  is  denoted 

by  ll/ll[-WW 

The  following  norm  is  described  in  more  detail  in  Section  2.2.  Let  jjZrt  be  the 
set  of  points  jVAf,...,iNAf-i-  This  notation  is  supposed  to  remind  the 

reader  of  the  set  of  points  {k/M}  modulo  N.  The  L2  norm  for  functions  defined 
on  jfZff  is  defined  to  be 


ll/lli>(wjv} =  E  \f{k/M)\7. 


k=-kJVM 


Our  version  of  the  FFT  is  an  operator  from  L2(—Zs)  onto  L3(^Zw)  defined 

by 

i  \NM~l  i  i 

Mn)  =  T7  E  f(k/M)e~i*iknfffM,  n  =  ijNTM-1. 

Mk =-\NM  2  2 

The  FFT  of  /  will  be  denoted  by  /. 

Note:  The  nonstandard  scaling  in  the  definitions  of  the  FFT  and  the  discrete 
norm  are  used  to  make  the  FFT  is  a  unitary  operator.  With  these  definitions, 

|l/|b(N,W)  =  I  I/I  b( Af.fr). 
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(See  Section  2.2  for  proof.) 

The  FFT  can  easily  be  extended'to  all  of  R  so  that  it  is  a  periodic  function 
with  period  M.  The  extended  FFT,  denoted  by  /*,  is 

.  yVAf-1 
fh)  =  77  E 

Note:  The  extended  FFT  is  also  a  unitary  operator.  By  Parseval’s  equality 
for  periodic  functions, 

il/B||[-Af/2,M/3]  =  ||/||d(W,W). 


The  periodization  operator  Py  is  defined  by 

/(^  +  &N),  /  6  L*(R)  and  decays  sufficiently  fast. 

kez 

See  Section  2.2  for  a  more  precise  definition. 
x  means  the  the  complex  conjugate  of  z. 

[zj  means  the  greatest  integer  less  than  or  equal  to  z. 


0.3  Executive  Summary 

FFT  Approximation 

From  now  on,  throughout  this  paper,  all  functions  are  assumed  to  be  in  L!(R), 
unless  stated  otherwise.  We  also  assume  that  /  and  /  are  continuous.  (In  fact, 
most  of  our  results  depend  on  the  existence  of  time  and  frequency  energy  moments, 
which  is  a  stronger  condition  than  the  continuity  of  /  and  /.) 

Because  of  our  assumptions,  both  the  Fourier  transform  and  the  FFT  exist. 
We  would  like  to  measure  the  error  in  approximating  the  Fourier  transform  by 
the  FFT. 

The  error  in  the  FFT  approximation  can  be  measured  in  many  ways.  One  way 
is  to  measure  the  maximum  difference  over  all  of  the  points  upon  which  the  FFT 
is  defined. 

Since  /  is  assumed  continuous,  we  can  evaluate  /  at  Therefore,  we  can 

also  measure  the  error  by  calculating  the  discrete  I ?  distance,  || /  - 
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We  can  also  measure  the  error  by  calculating  the  continuous  V  distance,  ||/«  - 

With  the  Ll  hypothesis,  a  key  step  in  deriving  an  error  bound  was  to  show 
that  /  and  /  decay  at  a  sufficiently  fast  rate.  We  obtain  a  similar  result  easily  in 
Lemma  1.5.  If  Cj  and  D\  are  finite,  then 

\f\t)\  <  A\t\-\ 

where  A  is  a  constant.  Unfortunately,  this  lemma  doesn’t  give  us  any  information 
about  what  the  constant  is. 

However,  a  different  approach  yields  an  even  stronger  result.  Define  the  dis¬ 
crete  tail  energy  to  be 

z  mm?- 

|Jfc|>lAW 

In  Lemma  2.3,  we  show  that  if  Cj  and  D\  are  finite,  then  the  discrete  tail  energy 
is  0{N  ■*),  and  we  also  have  an  upper  bound  for  the  constant  involved. 

Notice  that  this  implies  that  /3(t)  is  0(|l|--»).  At  this  point,  we  could  plug  this 
information  into  the  rearrangement  of  the  Poisson  summation  formula  as  we  did 
with  the  L1  hypothesis  to  obtain  a  bound  for  the  maximum  error.  This  approach 
results  in  the  first  FFT  approximation  theorem  (Theorem  5.1). 

By  preceding  somewhat  differently,  we  can  obtain  bounds  for  the  discrete  and 
continuous  L3  error.  The  following  briefly  describes  how  we  derive  an  upper  bound 
for  the  discrete  L3  error. 

Since  the  discrete  energy  in  the  tail  is  bounded,  one  would  expect  that  we 
can  find  an  upper  bound  for  the  discrete  V  distance  between  a  function  and  its 
periodization.  In  fact,  we  prove  the  following  periodization  comparison  lemma 
(Lemma  2.8).  If  Cj  and  Di  are  finite,  where  j  >  2,  then 

\\P#f  -  !\\d(M,N)  £  KN~J/2y 

where  K  depends  only  on  Cj  and  D\  (and  inversely  on  N  and  M). 

This  lemma  allows  us  to  replace  /  and  /  in  the  definition  of  the  discrete  L3 
distance  between  the  FFT  and  the  Fourier  transform  with  periodizations  (using 
the  triangle  inequality).  Of  course,  replacing  /  with  {PNfY  introduces  an  error 
of  KXN  :  and  replacing  /  with  P^f  introduces  an  error  of 

Now  we  must  calculate  \\(PNf)~  -  P„f\\D(NMy  For  this,  we  use  the  general¬ 
ized  Poisson  summation  formula  (Theorem  4.1).  The  generalized  Poisson  summa¬ 
tion  formula  says  that  under  the  hypotheses  of  the  Poisson  summation  formula, 
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(Prff)~  is  equal  to  the  sampled  Pm/.  Therefore,  in  our  case,  where  the  hypothe- 
■  ses  are  met  because  of  energy  moment  conditions,  the  norm  of  the  difference  is 
zero. 

We  have  just  derived  the  desired  upper  bound  for  the  discrete  L2  error  (The¬ 
orem  5.2).  If  Cj  and  Dj  are  finite,  where  j  >  2,  then 

\\i-f\\D(NM)  <  +  KiM-i/’. 

Ki  depends  only  on  Cj  and  D\  (and  inversely  on  N  and  M).  K*  depends  only  on 
Dj  and  C\  (and  inversely  on  N  and  M). 

We  also  obtain  an  upper  bound  for  the  continuous  L3  error  (Theorem  5.3,  see 
Section  5.3  for  the  proof).  If  Cj  and  Dj  are  finite,  where  j  >  2,  then 

II/'  -  flk-M/iM/2)  <  KiN-tl'  +  K,M-K 

K\  depends  only  on  Cj  and  D\  (and  inversely  on  N  and  M).  ivj  depends  only 
on  Dj.  Notice  that  the  second  error  term  is  better  than  the  corresponding  term 
in  the  previous  result. 


Riemann  Sum  Approximation 

Since  we  are  focusing  on  the  L3  norm  error,  our  approach  does  not  automatically 
provide  an  error  bound  for  more  general  Riemann  sum  approximation.  However, 
from  the  third  FFT  approximation  theorem,  we  derive  the  L3  norm  approximation 
theorem  (Theorem  6.1).  Then,  we  bootstrap  our  way  to  derive  error  bounds  for 
more  general  Riemann  sum  approximation. 

From  the  third  FFT  approximation  theorem  and  the  triangle  inequality,  we 
have 

|  ll/'llt-M/WJ]  -  II/II[-wam/>i|  <  KiN-H*  + 

As  mentioned  above,  the  first  term  of  this  inequality,  ||/*||(-\f/2,Af/2],  is  equal 
to  ||/[|o(Ar^r)  by  Parseval’s  equality  for  periodic  functions. 

The  second  term,  j|/||[_jv/2,Af/2),  can  be  replaced  by  ||/||  =  ||/||  at  the  expense 
of  increasing  K7.  This  is  a  consequence  of  Chebyshev’s  inequality. 

Making  these  two  substitutions  proves  the  L2  norm  approximation  theorem. 
If  Cj  and  Dj  are  finite,  where  j  >  2,  then 

|  ll/ll  -  [|/llz>(«j»)|  <  *1  N-«*  + 
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